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Chapter 1 



Introduction 

The outcome of this thesis is an algorithm solving the radial Dirac equation (to- 
gether with the corresponding computer code) in a specific case, as a part of a 
particular method for electronic structure calculations in solid-state physics. 

The electronic structure calculations are essential for any theoretical study of 
materials, which is currently an extremely active field of research and is often de- 
noted as computational material science. In past few decades, electronic structure 
calculations made a significant contribution to our understanding of material prop- 
erties, using a large variety of continuously developing methods and approaches. 

This work is related to one of up-to-date ab-initio pseudopotential methods 
based on the density-functional theory (DFT)[15, 9], particularly to the pseudopo- 
tential generating process within the all-electron pseudopotential method [13, 14]. 

Ab-initio method means that the method do not require any empirical param- 
eters as an input to the calculation — being based on equations derived directly 
from theoretical principles, with no need for experimental input data. 

The standard result of the particular calculation within DFT is the total energy 
of the given system and the electronic charge density (or the wave function of elec- 
tronic states that, integrated over the space and occupied states, form the charge 
density). Within DFT, the electronic charge density is a key quantity containing 
complete information about the system. This information allows us to answer al- 
most any question we might ask about the solid: the structure of electronic states 
provides information about thermal and electrical conductivity, e.g. metallic or 
semiconductive behavior, optical or X-ray emission/absorption spectra and scat- 
tering properties; minimization of the total energy with respect to atomic positions 
provides an equilibrium geometry; derivatives of the total energy gives bulk mod- 
ulus and elastic constants; spin structure of electronic states and derivatives of the 
total energy with respect to an external fields can provide dielectric or magnetic 
properties of the solid, etc. 

At present, several more-or-less ab-initio computational methods withing DFT 
are widely used in solid state physics. Every method usually has it's own domain in 
which it excels. The success of a method is determined by many factors, including 
its usage in particular computer codes and their efficiency, accuracy etc. And new 
methods can emerge in the future. It is not the aim of this thesis to compare all the 
methods available. It will suffice to say that in present, the majority of calculations 
are based on the local-density approximation (or its extensions) to DFT, which we 
review very briefiy now [9]. It leads to the following Kohn-Sham equations (for N 
electrons) : 

1 ^ 1 

_V2 + J2 K°"(x - X,) + Fh(x) + Fxc(x) \ M^) = e,^,(x) , (1) 



5 



where the charge density is 

N 

1=1 

the Hartree potential Vh is given by 

V^Vff = -47rp , 



and the exchange-correlation potential Vxc is just a given function of the charge 
density p. The is the electrostatic potential of atomic nuclei. Within the family 
of up-to-date pseudopotential methods, it is constructed from nonlocal pseudopo- 
tentials. There are several ways how to do that. 

These equations need to be solved self-consistently, so that the charge density, 
which is used for the construction of the Hartree and exchange- correlation potential 
is the same as the charge density calculated from the equations, see also fig. 1. 



initial p, 



V^V^ [p\^V^.[p] + Vi 



Hi) = Eip 

i ® 



converged to self consistency? 



yes 

Fig. 1: Self-consistency cycle 

Numerous methods have been developed to solve the resulting system of one- 
particle Kohn-Sham equations (1). Depending on the basis set used and other 
features of the particular method, some of them (in fact: almost all of them) need 
the radial part of the equation to be solved. For heavy atoms (approximately, 
starting at atomic numbers around 40-50) the equations (1) are considerably inac- 
curate due to the non-negligible relativistic effects, especially for the core electronic 
states bounded with relatively high energies. For this reason, it is desirable to use 
a relativistic Dirac equation instead, which is the aim of this thesis. 

In the present thesis we first derive and show how to solve the radial Schrodinger| 
equation for a given energy. Then we move to the relativistic case, Dirac equation. 
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and do the same. Next we show how to solve the eigenproblem, that is, how to de- 
termine the energies, for which the solution is normalizable. Then we say something 
about the computer implementation together with some results of the program. In 
the appendix we define atomic units (which is sometimes a little confusing issue) 
and we derive spin-angular functions including some of their properties we need in 
our work. 
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Chapter 2 



Schrodinger Equation 
2.1 Introduction 

The Schrodinger equation describes a nonrelativistic particle in a potential field. It 
cannot be derived, we always have to postulate something, more or less equivalent 
to the equation itself: 



where %l){x,y,z) is a (complex) function, p = —ihV, M is a mass of the particle 
(in our case an electron), V{x,y,z) the potential field (in our case we only have a 
spherically symmetric field V{r)). A physical quantity we are actually interested 
in is a probability density p = 

The electron has to be somewhere in the universe, thus we want 



For a given energy E, we can always solve the equation, but for the energies not 
lying in the spectrum of H, the solution exponentially diverges to infinity and such 
a solution cannot be normalized so that (3) holds. To be precise - there actually 
exist physical solutions, which are not normalizable according to (3) (the ones lying 
in the continuous part of the spectrum, for example a free electron, V = 0), but in 
our case of bounded states in a potential, we always have a discrete spectrum. 

The condition (3) picks up only certain energies (eigenvalues), when the solu- 
tion doesn't diverge. We label the energies by an integer number n starting from 
the lowest one n — 1, second lowest n — 2 etc. 

Besides energy, the solution also depends on quadrate of angular momentum (l) 
and it's z component (m). As it turns out, the radial part of the solution depends 
on n and I only. 

So we want to solve the eigenvalue problem of finding a solution for a given n 
and /. 

2.2 Radial Schrodinger equation 

We have a spherically symmetric potential energy 



State with a given square of an angular momentum (eigenvalue /(/ + !)) and its z 
component (eigenvalue m) is described by the wave function 




(2) 




(3) 



V{x) = V{r) . 




(4) 
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where Rni{r) obeys the equation [4] (eq. 2.400) 

Kl + -Kl + - y)Rnl - ^-^^Rnl = . (5) 

This is called the radial Schrodinger equation which we want to solve numerically. 

The derivation is well-known [4, 11], so just briefly. Basically, it's just a sep- 
aration of variables: we decompose the space as a tenzor product IR^ = x R, 
where is a unit sphere and IR is the radial part. We choose a basis in S^, it 
turns out that spherical harmonics Yim are a good choice as they are eigenvectors 
of and L3. We will search for all solutions of the form (4). Substituting (4) 
into the equation (2) yields (5): the trick is to write in spherical coordinates, 
the angular derivatives will then act on Yim only, thus separating the equaion. It 
turns out, that all the solutions Rni form a basis of R. So we have found a basis 
of R^, which is also a solution of (2) and thus any other solution can be found as 
a (possibly infinite) linear combination of V'n/m- 

2.3 Numerical integration for a given E 

Equation (5) is the linear ordinary differential equation of the second order, so the 
general solution is a linear combination of two independent solutions. Normally, the 
2 constants are determined from initial and/or boundary conditions. In our case, 
however, we don't have any other condition besides being interested in solutions 
that we can integrate on the interval (0, 00) (and which are normalizable) , more 
exactly we want i? e and Jq°° r^i?^ dr = 1. 

It can be easily shown by a direct substitution, that there are only two asymp- 
totic behaviors near the origin: and r^'"^. We are interested in quadratic in- 
tegrable solutions only, so we are left with and only one integration constant, 
which we calculate from a normalization. This determines the solution uniquely. 

All the integration algorithms needs to evaluate R", which is a problem at 
the origin, where all the terms in the equation are infinite, although their sum 
is finite. We thus start to integrate the equation at some small tq (for example 
ro = 10"^*^ a.u.), where all the terms in the equation are finite. If we find the initial 
conditions R{ro) and R'{ro), the solution is then fully determined. 

If ro is sufficiently small, we can set -R(ro) = Tq and R'{ro) = lr^(^^. This 
works fine for / > 0. For / = 0, it is not strictly correct, but it works well in practice 
because the fourth-order Runge-Kutta method is able to quickly correct the initial 
derivative guess. 

So when somebody gives us / and E, we are now able to compute the solution 
but the multiplicative constant that is later determined from a normalization. As 
was already mentioned, we used the fourth-order Runge-Kutta method that proved 
very suitable for this problem. 
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2.4 Asymptotic behavior 

The asymptotic behavior is important for the integration routine to find the correct 
solution for a given E. In this section we look into more details of the asymptotic 
expansion and illustrate it on 2 examples. 

It is well known, that the first term of the Taylor series of the solution is r\ 
independent of the potential [4] (eq. 2.408). This is enough information to find 
the correct solution for / > because the only thing we need to know is the value 
of the wave function and its derivative near the origin, which is effectively Tq and 
lr\^^ for some small tq. The problem is with Z = 0, where the derivative cannot be 
calculated just from I and tq. This section shows why and in the next section we 
show how we solved the problem. 

We start with the radial Schrodinger equation (5) and we shall search for the 
solution R in the form of a Taylor series: 

oo 

i? = ao + a\r + a2r^ + . . . = ^ akV^ . 
Substituting this into the equation we get: 

oo ^ oo 

^ T^au [k{k +1) -1(1 + 1)]+ {E-V)J2 ^'«fc-2 = . (6) 

k=0 k=2 

Let's assume we have a potential V of the form: 

oo 

V = \-vo + vir + V2r + . . . = ^ vjr-' , 

we rearrange the double sum on the right hand side of (6) 



oo oo oo oo 



k=2 j=-ifc=2 j=-ik=j+2 

oo oo oo fe oo k 

= YY1 Vj-ir^^^(^k-j = YY1 Vj-ir''+^(^k-j = Y = 

j=0 k=j k=0 j^O k=0 j^O 

oo k—1 

k=l j=0 



So we get: 



Y rW [k{k + 1) - l{l + 1)] + -^E Y r'afc_2 + 

k=0 k=2 
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2M 
If 



oo k — 1 



fe=l j=0 

This equation holds for every r, thus we coUect the coefficients at and they must 
vanish. We get: 



/c = aQ[-l{l + 1)] = Q , 



2M 



k = l ai[2- 1(1 + 1)] 2-v_iao = 0, 



(7) 
(8) 



k>2 ak[k{k+l)-l{l + l)]-^y"^Vj_iak-j-i + '^Eak-2 = 0. (9) 

This enables us to calculate all the Taylor coefficients of the solution. To 
see how it works, we calculate two examples and compare them to the analytical 
solution. First, let 

r 

so V-i = —Z, vq = vi = . . . = 0. For Z = 0, we see from (7) that ao can by any 
number (including zero, but as we will see in a moment, this would imply the zero 
solution, which we are obviously not interested in). From (8) it follows: 

MZ -ao 
ai — — r2~oo 



a 



where a = Bohr radius. The ffist two terms of the solution are then: 

i? = ao(l - ^ + 0(r2)) , 
which is in agreement with the analytic solution [4] (eq. 2.524) (every R^i for / = 0): 

-Rio 




2\l ^7 exp (-^ 



a 



R20 
R30 
R40 



1 



2a3 
2 / 1 



1 



2a 



exp 



2aJ ' 



3 V 3a3 
1 



^ 2r 2 / r \ 2 
~ 3a ^ 27 Va. 



1 



4a2/2 V a 




3r l/r\2 1 /r\3 

1-^ + 77 - - 



192 \a 



exp 



4a 8 Va 

As the second example, we use a linear harmonic oscillator 

' — 1 
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so v_i = vq = vi = 0, V2 = = = . . . = 0. For Z = 0, we see from (7) that 

ao is any number, from (8) it follows ai = and from (9) we get 02 = 
we know the spectrum [4] (eq. 2.484): 



MEap 



But 



E = nu;{2n + l + -) , 



so we have 02 



Mc<;(2n+|)ao 
3h 



— (|n+ 1)^5-, where we used the substitution 



a = 



\J~^^- Finally the first two nonzero terms of the solution are: 



R = ao[l- [^n + \] ^+0(r=^) 



which agrees with the analytic solution [4] (eq. 2.488) (again every R^i for Z = 0): 



00 



-Rio — 




exp 



/ r2 - 




2a^ 


-( 


r \ 2 


3 V 


a) 



cxp 



2a? 



These examples show, that for Z = the derivative R' (the second term in 
the R expansion) nontrivially depends on V in the first example, and on E in the 
second example. Which is inconvenient for a numerical computation. 

For Z > 0, the Taylor coefficients can be calculated in the same way as for 
I = 0. From (7), (8) and (9) we see that Ofc = for all k < I. So indeed the first 
nonzero term is air'' as expected. 
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Chapter 3 
Dirac Equation 

3.1 Introduction 

The Dirac equation for one particle is [12, 16]: 

HiP = Wtl;, (10) 

H = ca p + j3mc^ + V{r)t , 
where t/^ is a four component vector: 




and a, P are 4x4 matrices: 



P 



cr 

o- 

1 

-1 



where the Pauh matrices cr = {crx:(^y:'^z) and 1 form a basis of aU 2 x 2 Hermitian 
matrices. Substituting all of this into (10) yields: 

^^\^a)-[ W-V + mc')[i;B) ■ (11) 

To derive a continuity equation, we multiply (10) by i/j* and subtract the conjugate 
transpose of (10) multiplied by i/j: 

so we identify the probability and current densities as 

The normalization of a four-component wave function is then 

j pd^x = J V'^ = J V'l^i + ^p2^2 + V'sV'a + d^x=l. (12) 

The probability density p{x,y,z) is the physical quantity we are interested in, and 
all the four-component wavefunctions and other formalism is just a way of calculat- 
ing it. This p is also the thing we should compare with the Schrodinger equation. 
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3.2 Derivation of the radial equation 

We rewrite cr ■ p [2] : 

lcr-x/.9 ^\ 

(T • p = —inr— + zcr • L , (13) 

r r \ or J 

and search for a basis in the form 

= gxi' , (14) 

V'B^^/xi- (15) 

We chose this form, because we want i/^a and i/jb to have the same j and j3, so 
the only thing they can differ in is / (the two aUowed possibihties are I — j ± 
According to (36), (37), (38) and (37), this corresponds to ±k, (see also [6, 12, 10, 
2] and the appendix for more information about spin-angular functions). 
Substituting (14), (15) and (13) into (11) we get: 

r r \ Or J \ QXk J 

W -V -mc^ \ f gx^ 

W-V + mc^ J yifx-^ 

Both ipA and ipB are eigenvectors of the operator K = I3{cr ■ h + h): 

KiiA,B = -hKlljA,B , 

SO the action of the a ■ L operator on the i/ja and i/jb can easily be determined: 

..L(^^)^(,^-a)(;^^ 

K-n \f^^\^fn{-K-i)ti;A 
-K-hj {iJBj V M«-i)V'B 

and we can write 



c— 



1^-x / i-ihr^ + z(« - l)n)z/xi^ 
r r V(-^^^^ +^(-« - 1)%X:^' 



W -V -mc^ \ / ^xi' 

W-V + mc^ J \ifX-K 

The operator also only acts on the angular momentum parts of the state [12] 
(page 59). From (35): 

K A — K 5 

r 

14 



1 fithrl-ziK-mifxi' 



c 



dr 
A. 
dr 

W -V -mc^ \( gxi^ 

W-V + mc^ J U/xi' 



rewriting 



r\ (r-|: + («+l))^A^x'l 

^f{W-V-mc^)g \ / xi^ 

V iW-V + mc^)f J \ix'^^ 

and canceling the same terms on both sides finally yields 

V ^ + ^9 )-[iW-V + mc')f) ■ ^^^^ 

This is the radial Dirac equation. As we shall see in the next section, the equation 
for g is (with the exception of a few relativistic corrections) identical to the radial 
Schrodinger equation. And / vanishes in the limit c — > oo. For this reason / is 
called the small (fein, minor) component and g the large (groB, major) component. 
The probability density is 

p = ^> = rAi^A + rB^B = fx'r^x'x + 9\rx'^ , 

so from the normalization condition (12) we get 

j pd^x^ J fx'^^^x'^^ + g'xrxi'd'x^ J{fx'^:x'^. + 9\rxi')r'drdn^ 

= / fr^ dr / xi':X-« dO + / g^r^ dr / xi^*xi^ dn = / r^{f + g^) dr ^ 1 , 
Jo J Jo J Jo 

where we used the normalization of spin-angular functions (32). Also it can be 
seen, that the radial probability density is 



Pir)=r'{r + g') 

(i.e., the probability to find the electron between ri and r2 is J^^ r^{f^ + g^)dr). 
In the nonrelativistic case, the density is given by 

p(r) = r^R^ , 

so the correspondence between the Schrodinger and Dirac equation is E? = f'^+g'^. 
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3.3 Numerical integration for a given E 

We rewrite (16) to a form found in [12] (eq. 8.10 and 8.12): 

g'. = -^9. + ^{W-V + mc^)h , (17) 
r cn 

f. = "'^'^^'^ • ^^^^ 

Let h = 1 and define E = W — mc^, so that E doesn't contain the electron rest 
mass energy. Next we define a relativistic mass 

M = m + ^{E-V) (19) 

and we get 

^« = -^5'K + 2Mc/fe, (20) 

f'. = +'^U + -iV-E)g,, (21) 
r c 

which are the same equations as in [7] (eq. 2a and 2b). 

From (20) we express and substitute it into (21). At the same time we 
introduce a new variable ^ 

0. = , (22) 

(beware, [7] introduce 0^ = jWcd'^.^ simplify our equations: 
f = ^''^ K + 1 K + 1 

^'^ 2Mc r 2Mc c c 2Mr ' ^ ^ 

differentiate 

•'^ c c \2Mr) c cr 2M^cr^^2c'^ 2Mcr^^^ 

and substitute f^, and back to (21): 

K-1 f (1)^ K+1 g^ \ 1 



9.{r) . (24) 



r \ c c 2Mr ) c 
by a simplification we finally get 



^«(O = -^0K(r) + 
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Equations (24), (22) and (19) are the equations we are solving. It is instructive to 
write equation for g^^ directly. For this, we need to calculate 



9' 



g'M' 
2M ~ 2M2 

and substituting for and (f)'^ into (24) gives 



9'^ 



9'.M' 2 g'^ 



+ 



2M 2M2 r 2M 
multiplying by 2M and rewriting 

'2 M' 



iy -E) + 



K,{k + 1) K + 1 



2Mr2 4M2c2r 



9k 



M 



«;(«; +1) K + 1 
2Mr2 



4M2c2, 



Prom (19) it follows 



M' 

'm 



V 



2Mc2 



so we finally get 
'2 



^ 2Mc2 )^^^ 



k{k + 1) K + 1 



2Mr2 



4M2c2r 



2M^, . (25) 



Comparing (25) with (5) we clearly see the two relativistic corrections, both de- 
pending on V and both vanishing for c ^ oo as expected. Also it should be noted 
that there is another difference, that (25) contains the relativistic mass (19), but 
(5) contains the rest mass m. 

All the terms in the (25) are the same for both possibilities j = I ± ^ (i.e. 
k{k + 1) = 1(1 + 1) for both k = —l — l and k = I), except for the spin-orbit coupling 
term 4^T^V'{r), see for example [12] (eq. 2.64). Sometimes it makes sense to 
consider a semirelativistic case, where we neglect the spin-orbit coupling term, in 
which case we are left with the Schrodinger equation with the relativistic mass M 
and only one correction 2Mc^ ■ 

In practice, the potential V is given on a discrete grid, so we need to com- 
pute it's derivative V' numerically. This is the reason we solve (24) instead of 
(25), because we need to evaluate V only once for the spin-orbit term (for the 
semirelativistic case we don't even need the V at all) . But besides this minor tech- 
nical thing, there is no other reason we chose (24) and not (25), which looks more 
familiar. 

Once we have calculated g,^ and we calculate from (23). So the result 
of the radial Dirac equation are two functions and g^- The physically relevant 
quantity is the radial probability density 



Pir)=r\f,+9l) 
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We calculate the functions /«; and g^, in a similar way as we calculated R for 
the Schrodinger equation, thus we need the asymptotic behavior at the origin. The 
potential can always be treated as V — 1/r + ■ ■ ■ and in this case it can be shown 
[16], that the asymptotic is 

.13-1 



where 



{(3 - l)r^-2 



2M 




(26) 



or, if we write it explicitly, ioi j — I -\- \ 



P' 



+ 



v'h-i)^-( 



and j — I — 



2 



z 



2 



In the semirelativistic case (which is an approximation — we neglect the spin-orbit 
coupling term) we choose 



p = + i/3-p) = y p + z + 1 - (I) . 

It should be noted that in the literature we can find other types of aymptotic 
behaviour for the semirelativistic case, its just a question of the used approximation. 
One can hardly say that some of them are correct and another is not since the 
semirelativistic (sometimes denoted as scalar-relativstic) approximation itself is not 
correct, it's just an approximation. 

It follows from (26) that for j = Z + | the radial Dirac equation completely 
becomes the radial Schrodinger equation in the limit c — > oo (and gives exactly the 
same solutions). For j = l — \ however, we get a wrong asymptotic: we get a radial 
Schrodinger equation for Z, but the asymptotic for Z — 1. 
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3.4 Other forms of the equations 

Unfortunately, there are a fot of different forms the radial Dirac equation can be 
found in the literature. Many authors use different symbols, different units, I even 
found a mistake in [8]. For the reader's comfort, this section is devoted to deriving 
and presenting the most common forms of the equations. 
We realize the fact 



r \dr r 
, K + 1 1 / d k\ , , 



and rewrite equations (17) and (18): 



r \ar r J c 

-(^--) (^/-) + -{W-V- m^)gu = . 
r \dr r J c 

These equations could be found in [16] (eq. 8.10 and 8.9). Let's make the substi- 
tution [8] 

Qk — ^fn 

and write 

^ + -)p.--iW-V + mc^)Qk = 0, 
dr r / c 

^--]q. + -{W-V- mc^)Pk = , 
dr r J c 

which can be found in [3] (eq. 3 and 4: they write anijir) = and bnij{r) = Q^) 
and also in [12] (eq. 8.13: he uses = Pk and = Qk)- 
Now we use (19) and these relations become 

W -V + mc^ = E -V + 2mc^ = 2Mc^ , 
W-V-mc^^E-V, 



to write 



dP K 

-r^ = Pn + 2McQk , 

dr r 

^ = -Q.--{E-v)P,, 

dr r c 
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which can be found in [8] (eq. 2a and 2b) (there is a mistake there, they forgot to 
divide by r). 2M can be written exphcitly as 



2M 



E-V 




'E-V 




^ +2m 






a.u. 



so we get 



dr 
dr 



K 



= Pk + 



E-V 



+ 2 



^Q^-l(E-V)Pk 
r c 



cQh 



which can be found in [16] (eq. 8.12 and 8.13), where they have one c hidden in 
Qk = cr/«; and use Rydberg atomic units, so they have 1 instead of 2 in the square 
bracket. It can be found in [1] as well, they use Hartree atomic units, but have a 
different notation Gf^ = Pf^ and F,^ = Q^, also they made a substitution c—^. 
Some authors also use e = E. 
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Chapter 4 
Eigenproblem 



4.1 Introduction 

In the previous two chapters, we learned how to calculate the solution of both radial 
Schrodinger and Dirac equations for a given E. For most of the energies, however, 
the solution for r — > oo exponentially diverges to ±00. Only for the energies equal 
to eigenvalues, the solution tends exponentially to zero for r — > 00. The spectrum 
for bounded states is discrete, so we label the energies by n, starting from 1. 

We want to find the eigenvalue and eigenfunction for a given n and / (and a 
spin in the relativistic case). The algorithm is the same for both nonrelativistic and 
relativistic case and is based on two facts, first that the number of nodes (ie. the 
number of intersections with the x axis, not counting the one at the origin and in 
the infinity) of Rni and ^r^ is n — Z — 1 and second that the solution must tend to 
zero at infinity. 

4.2 Algorithm for solving the eigenproblem 

We calculate the solution for some (random) energy Eq, using the procedure de- 
scribed above. Then we count the number of nodes (for diverging solutions, we 
don't count the last one) and check, if the solution is approaching the zero from 
top or bottom in the infinity. From the number of nodes and the direction it is 
approaching the zero it can be determined whether the energy Eq is below or above 
the eigenvalue E belonging to a given n and /. The rest is simple, we find two 
energies, one below E, one above E and by halving the interval we calculate E 
with any precision we want. 

There are a few technical numerical problems that are unimportant from the 
theoretical point of view, but that need to be solved if one attempts to actually 
implement this algorithm. One of them is that when we end the algorithm, because 
the energy interval is sufficiently small, it doesn't mean the solution is near zero 
for the biggest r we are integrating the equation. Remember, the solution goes 
exponentially to ±00 for every E except the eigenvalues and because we never find 
the exact eigenvalue, the solution will (at some point) diverge from zero. 

Possible solution that we have employed is as follows: when the algorithm ends 
we find the last minimum (which is always near zero) and trim the solution behind 
it (set it to zero). 

Another solution of this problem is to end the algorithm not only when the 
energy interval is small enough, but also when the solution is sufficiently near zero 
for the biggest r we integrate. One would thought (including me at first :) that 
in this case we don't need to trim the solution. Mistake. If the biggest r is big 
enough, then even the variation of the initial E by 10~^^ is not enough to push 
the tail of the solution to zero. It can actually be above zero by even 50% of the 
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maximum of the wave function or more, so we would have to trim the solution 
anyway (not mentioning that the algorithm never ends, because the solution will 
never be close enough to zero to pass the ending condition — that actually cannot 
even be formulated. . . ). 

The second rather technical problem is how to choose the initial interval of 
energies so that the eigenvalue lies inside the interval. 
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Chapter 5 
Implementation 



5.1 Fortran routine 

The original Fortran 77 routine solves the equation (5) on a hyperbolic grid (27) 
using a simple polynomial approximation at each step. The input is the energy E, 
grid parameters ap, ko and Jm, plus an array with values of the potential V{r) at 
the grid points. The output is an array with values of REi{r) at the grid points. Our 
task was to replace this routine (with minimum changes in other parts of the code) 
with the new routine that solves the equations (24), (22) and (19) on the hyperbolic 
grid using the fourth-order Runge-Kutta method. The input is the same as for the 
original routine, plus Z. The output are the two functions / and g. 

This means that in the limit c — > oo the g component of the new routine gives 
the same results as the original one. 

This new routine can switch between the fully relativistic (Dirac) case, a 
semirelativistic case (we neglect the spin-orbit coupling term) and a nonrelativistic 
case (we neglect the spin-orbit coupling term and use a nonrelativistic mass, which 
effectively restores the Schrodinger equation). 

5.2 Eigenvalue problem 

This routine uses the algorithm described in the preceding chapter. It can use 
either the new or the old routine for the integration. 

5.3 Hyperbolic grid 

The routines solve the equations on a hyperbolic grid, which is defined as 

^ = «^l4^' (27) 

where < ^ < 1 is a dimensionless parameter and ap is a scaling (in the same units 
as r), it is the value of r for ^ = 0.5. We take equidistant values for ^ 

C= ko<jM, (28) 

jM 3m Jm 

by substituting (28) into (27) we get 

r{j) = ap-^, j = l,2,...,A;o. (29) 

JM-J 

The following values of the grid parameters has proven successful: ap = 1, ko = 500 
and Jm = 525. 
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Routines accepts the parameters ap, ko and Jm, then the array r{j) of a 
dimension ko, which gives the values of r in the points j = 1, 2, . . . , /cq calculated 
using the relation (29). The routine can easily calculate all the values of r(j) from 
(29), or use the values directly from the array r{j). We use the way that is more 
readable. 

A formula for the inverse transformation can be derived from (29): 

3 = 3m- 
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Chapter 6 
Results 



6.1 One electron atom calculation 

To check that our program works, we have computed a spectrum of a one electron 
atom with these parameters: Z = 92, ko = 2920, Jm = 3000, ap = 1.0. The result 
can be compared to fig. 8.3 in [12]. 

In the following table, Ei represents eigenvalues calculated from the nonrela- 
tivistic formula 

2n2' 

E2 represents eigenvalues in the nonrelativistic case (so it should not depend on /) 
and E3 eigenvalues computed from the radial Dirac equation, so it depends on I. 



The energies are given in eV. 



Case j — 


i + h 




















n=l, 


1=0, 


El=- 


■115156, 


.9520 


E2=- 


-114712, 


.4604 


E3=- 


-131994, 


,5289 


n=2. 


1=0, 


El= 


-28789 


.2380 


E2= 


-28733, 


.2317 


E3= 


-34151, 


,9041 


n=2. 


1=1, 


El= 


-28789 


.2380 


E2= 


-28789, 


.2246 


E3= 


-29649, 


.2912 


n=3. 


1=0, 


El= 


-12795, 


.2168 


E2= 


-12778, 


.5942 


E3= 


-14649, 


.7703 


n=3. 


1=1, 


El= 


-12795, 


.2168 


E2= 


-12795, 


.2121 


E3= 


-13307, 


.1396 


n=3, 


1=2, 


El= 


-12795, 


.2168 


E2= 


-12795, 


.2167 


E3= 


-12959, 


.5541 


n=4. 


1=0, 


El= 


-7197 


.3095 


E2= 


-7190, 


.2921 


E3= 


-8026, 


.0999 


n=4, 


1=1, 


El= 


-7197 


.3095 


E2= 


-7197, 


.3073 


E3= 


-7466, 


.8873 


n=4, 


1=2, 


El= 


-7197 


.3095 


E2= 


-7197, 


.3094 


E3= 


-7318, 


.8303 


n=4, 


1=3, 


El= 


-7197 


.3095 


E2= 


-7197, 


.3094 


E3= 


-7248 , 


.7232 


Case j = 


/ - i (for / 


= we take j 


= h) 












n=l, 


1=0, 


El=- 


-115156, 


.9520 


E2=- 


-114712, 


.4604 


E3=- 


-131994, 


,5289 


n=2. 


1=0, 


El= 


-28789, 


.2380 


E2= 


-28733, 


.2317 


E3= 


-34151, 


,9041 


n=2, 


1=1, 


El= 


-28789 , 


.2380 


E2= 


-28789 , 


.2246 


E3= 


-34226, 


.2202 


n=3. 


1=0, 


El= 


-12795, 


.2168 


E2= 


-12778, 


.5942 


E3= 


-14649, 


,7703 


n=3, 


1=1, 


El= 


-12795, 


.2168 


E2= 


-12795, 


.2121 


E3= 


-14673, 


,2066 


n=3, 


1=2, 


El= 


-12795, 


.2168 


E2= 


-12795, 


.2167 


E3= 


-13307, 


,1884 


n=4. 


1=0, 


El= 


-7197 


.3095 


E2= 


-7190, 


.2921 


E3= 


-8026. 


,0999 


n=4, 


1=1, 


El= 


-7197 


.3095 


E2= 


-7197, 


.3073 


E3= 


-8035 , 


.9706 


n=4. 


1=2, 


El= 


-7197 


.3095 


E2= 


-7197, 


,3094 


E3= 


-7466 , 


.9097 


n=4, 


1=3, 


El= 


-7197, 


.3095 


E2= 


-7197, 


.3094 


E3= 


-7318, 


.8305 



It should be noted that the eigenvalues don't depend on I in the nonrelativistic 
case only for exactly the Coulomb potential. In practice, for potentials screened by 
other electrons in real atoms, this is not the case. 
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Chapter 7 



Conclusion 

Radial numeric solution of the Dirac equation has been implemented under specific 
conditions for all-electron pseudopotential genereitmg process and the corresponding 
computer code has been debugged. 
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Chapter 8 
Appendix 



8.1 Atomic units (au) 

There are two types of atomic units — Hartree units and Rydberg units. The term 
"atomic units" usually means Hartree atomic units and we use it in this meaning 
as well. Some authors however [16] mean Rydberg units when they write "atomic 
units" . . . 

All (Hartree) atomic units can be derived from this relation: 

h = m = e = Atteq = 1 , 

where m is the mass of the electron, e a charge of the electron, eo the permittivity 
of vacuum. 

Examples: dimension of length is a Bohr radius 

ao = = J!- = 1 a.u. = 0.529 A = 0.529 • lO-^o m , 

where a is the dimensionless fine-structure constant 

e2 1 

a 



ATTeohc 137.036 ' 

so the speed of light is 



c—— a.u. = 137.036 a.u. . 
a 



Energy is measured in Hartrees, one Hartree is 



Eh = n = mc^a^ = 1 a.u. = 27.211 eV . 



Hydrogen atom Hamiltonian in SI units: 

. 1 



H = 



2m 47reo r 
and au units: 

2 r 

The energy spectrum for hydrogen is 



E. 



moQ 2-n?' 2n2 
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The Rydberg atomic units are defined as 



in these units 



h = 47reo = 1, m = - , e^ = 2. 



ao = 1 , 

2 2 1 

Eh = 2 = "^c a = - , 

moQ 2 

r 

1 1 



E. 



c = - = 274.072 . 

a 



This is why the asymptotic term Za is mostly written as ^ (Hartree units), but 
sometimes as ^ (Rydberg units) . 

8.2 Spin-angular functions 

The operators Lp' (quadrate of angular momentum), L3 [z component of angular 
momentum) and {z component of a spin, ±|) mutually commute, so let's denote 
the corresponding normalized eigenvectors by \lmss). We can however also use 
commuting operators (quadrate of total momentum), and J3 {z component 
of total momentum) and denote their normalized eigenvectors by Ijljs)- 

First we want to find a relation between these two bases. A standard procedure 
is to use Clebsch-Gordan coefficients [5] : 



\jij3) = ^ X] h ™' *3 1 J, is) \ims3) = (30) 

m=-l s3 = -5 

= {I, |, js -hi 1^'' ^3) K;j3 - i; i) + h ^3 + |, li, is) \i;j3 + h -1) ' 

where (ji, j2, mi, m2 | J, m) are Clebsch-Gordan coefficients and we used the fact 
that they are nonzero only for rrii + m2 = m. Because wc arc adding angular 
momenta / and |, the only possible values of j are j = / ± i. j = I + k means, that 
the spin has the same orientation as the angular momentum, j = I — corresponds 
to the opposite orientation. So we will only need these 4 coefficients (see eq. K58 
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in [4] or table 2.2 in [12]) 

21 + 1 



{I, ij3 + i-|U + ij3) = 

(^) 2' -^3 ~ 2' 2 I ^ ~ 2' -^3) = 



2 ' -^'3 + 2 ' 2 I ^ 2 ' -^'3) 



21 + 1 ' 

'^ -J3 + | 

2/ + 1 ' 



2Z + 1 

Substituting into (30) yields the desired relation between the two basis: 



li = ^ - I; hh) = -\l ^^^^ ^ ^ = i3 - |; |) + 



Now we want to find representations of these "kets". Our Hilbert space can 
be decomposed as 

= ^ C^ 

where is a vector space of quadratic intcgrable functions on a unit sphere (the 
angular momentum operator L acts here) and is a two dimensional complex 
vector space (the spin operator 5" acts here). 

The eigenvectors of the S3 operator form a basis of the space C^, let's denote 
them where s = i^: 

The eigenvectors of and L3 form a basis of S^, they are called spherical harmonics 

Ylm. 
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Then 



\lm;s3 = I) = Yim^i = 
\lm;s3 = =Yim^_i ^ (^y^ ^ 



\j = i- hih) = -^^^ (-Vi-i3 + |y,,^=,3-i^i+ 



2 1 _ „ J3 

- Vji 



The yj^ are caUed "spin-angular functions" [11], or "spin spherical harmonics" [16]. 
For the two possibilities j = Z ± | we introduce a new quantity k: 



K 



/ — 1, for j = Z + i; 
Z, forj = Z-2, 



so that we can deduce from k both / and the spin orientation and thus we can label 
all the spin-angular functions by just two parameters n and j3 and write: 

Any integer value of k is permissible except k = 0. Some authors also use a different 
notation js = rrij = n for the eigenvalue of J3. 

It can be easily shown, that because of the normalization of spherical harmonics 

Yim- 

^ YC^{p)Yi,m'{p)d^ = 5u'5r, 
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the spin-angular functions have the following normalization: 

/ X{^*(n)x5(n)dl7 = 5««,5,3,.. (32) 

To be more precise, —Hk is actually defined as the eigenvalue of the operator 
[12] (eq. 2.85 and 2.86) 



K^= |^J^(j2_L2-S2) + l^nV' = 

= (j(j + 1) - /(/ + 1) - s{s + 1) + 1) = 

- {j{j + i)-i{i + i) + \)n^ = 

= —Hki/j , 

so we have k = + 1) + l{l + 1) — \. For j = I + ^: 

K = -{I + + ^ + 1) + + 1) - I = -I - 1 
and for the other case j = I — \ : 

K = -{l-\){l-\ + l) + l{l + l)-\ = l. 

Another way of looking at k follows from 

r_;_i = _(_^-+ 1), forj = / + i; 
U = +(J + |), forj = Z-|, 

which means that \k\\s actually the value of the total angular momentum (plus ^ 
so that we only have integer values) and the sign depends on whether the spin and 
the orbital angular momentum have the same or opposite orientation (j = Z ± |). 
In the text, we need to find how the operator 



acts on the spin-angular functions. It commutes with J, as is easy to prove, which 
means it also commutes with J3 and and so ^xtt must have the same values 
of j and j3 as x^k- Furthermore, it anticommutes with K [12], so 

K^xi^ = -^Kxi^ = -^{-nK)xi^ = hK^xi^ , 

r r r r 
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but also 



which shows that 



(33) 



for some (complex) constant C. From now on, there are mistakes in the literature. 
For example [12] implicitly concludes at this point that C = — 1, which is obviously 
incorrect. It can easily be shown that (^^)^ = 1, so = 1, or C = ±1. The sign 
of C can generally depend on k and j3, but doesn't depend on x. To determine it, 
one option is to find it explicitly for all values of k and J3. [10] chooses x = e^, the 
unit vector along the 2;-axis, and writes 



Yim{ez) 



21 + 1 
47r 



(34) 



<T • X 



SO using (33) 



3,3±2 3,JT2 



where the + sign corresponds to k > and the — sign to k < 0. Using (31) 

U ± i i J3 - i +k \j, h)Y 1 1 

-u ± i i is + i -| I j, j3)y.^i^^^^ 

U T i i is - i j3)Y.^i^.^_i 



(i^i ii3 + |,-||i,is)i;.^i .^^1 



evaluating the Clebsch-Gordan coefficients 



1 



V2i + 1 ± 1 



1 



c 



V 



3To,J3- 



1 . , 1 / 

T2.J3+2 / 



Using (34) it follows that C = — 1 for j > and J3 = ±i (and both ±). For other 
values of js however, the spherical harmonic Yiyn{ez) — (because m ^ 0), so we 
can say nothing. Unfortunately, js can be any number from —j, —j + l,...,j,so 
proving that C = — 1 only for js = ±| is insufficient and [10] is wrong as well. 
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We'll give a correct proof using a Wigner-Eckart theorem that C = — 1 indeed in a 
moment. So we finally get 

■ = -Xi^ • (35) 



r 

J3 



Let's look at the relation between and X-k niore closely. For k>0: 

Xi^ = ei,.' (36) 



for K < 0: 



Xk'' — 1 , — , ' (38) 

Basically, x^k X-k have the same j and js, and differ in I. There are only 
two possibilities I = j ± ^, so X'li picks one and X-k other one (which one 

depends on the sign of k) . So the operator preserves j and js and fiips the I to 
the other possibility. 

Let's get back to determining the constant C. We denote the eigenvector of 
Ji , JI , and Jz by a ket 

\j1j2jm) , 
where J = Ji + J2. In this notation, 

Xk =(x|j±^;|jm) , 

where Ji is the orbital angular momentum, J2 is the spin, j — \n\ — \ and the 
plus sign is taken for k > 0, negative sign is for k < 0. We want to determine the 
constant C in 



that is 



a ■ X 



r 

The kets \jij23T^) are orthogonal, so 



^x{^ = Cx'l^ , 



j±l;ljm) = C\jTl;ljm) 



c = 0' T |; — li ± ^; li"^) • 



We introduce two vector operators 



r 
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and form tensor operators T(l, g) and q) {q = —1, 0, 1) in the standard way: 

r(i,o) = ^, T(i,±i) = T^^, 

V(1,0) = ,T3, V(l, ±1) = T^i^^ , 

they obey the foUowing commutation relations 

[T(l,g),J2] =0, 
[V(l,g),Ji] =0, 

so we can use the Wigner-Eckart (W-E) theorem for the scalar product of two tensor 
operators (note that T(l)-V(l) = ^^{-1)1T{1, q)V{l, -q) = TiVi+T2V2+TsVs): 

C^ijT h iM^lj ± i; ijm) = {j T h ^^'^1^(1) ■ ± i; ijm) = 

= (_l).+.±i + ||i i|i _^.|i}(jT|||T(l)||j±i)(i||V(l)||i). 
Now (remember j can be either integer or half integer) 



1 ^ .li = . ^ 1 i ^ 



1 2 i±2j U 2 i ^6(j + \) 

(jT|l|Tr(i)||j±|) = TVJ+i, 

(il|V(l) ||i) = v^, 



so 



The 6j symbol was determined from the formula [4] (K.108) 



a b c\ _ ^_^^a+b+c, li(^ + b + c+l)i-a + b + c)(a-b + c) 



1 c-1 bj ^ ' y 6(26+ l)(26 + 2)(2c- l)2c(2c+ 1) 
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■y/a + b — c+1 , 

The reduced matrix element ||V(1) \\^) was determined from another apphcation 
of the W-E theorem: 



rnvMiu) = (-i)^^-i^i^MMa^(i II v(i) III) 

'2-k + l 



but also 



= ^(|l|V(l)||i), 



so 



(i||V(l)||i) = V6. 
Similarly for the reduced matrix element (j ^ | || T(l) || j ± |): 

{j T |;m|T(l,0)U± i;m) = ^-iy^k+^-^Tk0^i±i^!^^il^ x 



X (i T I II T(i) II i ± i) = ^ ^ II II ^' ± ^) ' 

but also 

{j T i;m|T(l,0)|j ± i;m) = (j ^ ± i;™) = 



1 /('o■ + i^2_™2 

y* 1 (i?,^)cosi?r.^i. (i?,^)df]=-^/^-^ 2^ 



so 



(jT|||Tr(l) ||j±i) = T^/J7i. 
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